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A  Large  Eddy  Smulation  (LES)  was  used  to  model  a  two  phase  turbulent  round  jet.  The  LES  for 
e  continuous  phase  used  the  Smagorinsky  model  for  the  sub-grid-scale  motion.  A  simple  model 
tor  me  effect  of  me  particles  on  the  continuous  phase  was  developed.  The  model  was  essentially  a 
momenrnm  balance  for  a  fluid  volume  containing  particles  moving  relative  to  the  continuous 
pnase.  ihe  simulation  showed  that  particles  accelerated  the  roll-up  of  the  shear  layer.  Droplet 
dispersion  was  measured  in  a  turbulent  round  jet  of  air  issuing  into  still  air.  Individual  droplets 
were  formed  with  a  piezo-electric  device;  they  were  injected  onto  the  centerline  of  the  iet.  Their 
locauon  across  the  jet  was  measured  at  different  axial  stations  with  a  laser  sheet  and  position 
sensing  photomultiplier  mbe.  Particles  containing  a  fluorescent  dye  were  injected  into  a  turbulent 
spray,  ^pt^  scattering  was  rejected  with  a  Raman  filter.  Dispersion  of  the  tagged  particles  was 
me^ured.  Results  agreed  with  the  LES  in  which  the  injection  of  the  spray  enhanced  development 
of  me  jet  and  me  dispersion  of  particles. 
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1.0  Executive  Summary. 


Experiment. 

Particle  dispersion  and  particle  velocities  were  measrued  with  laser  sheets  and  a  position 
sensitive  photo  multiplier  tube  to  track  particles.  Monodisperse  hexadecane  droplets  were 
injected  onto  the  centerUne  of  a  tiubulent  air  jet.  Their  radial  dispersion,  axial  velocities, 
and  times  of  flight  were  measured  as  functions  of  axial  position.  The  time  and  length 
scales  of  the  jet  were  varied  through  the  control  of  the  jet  exit  velocity  and  nozzle  diam¬ 
eter.  Nozzle  diameters  of  7  fim  and  12.6  fxm  were  used.  Reynolds  numbers  were  in  the 
range  of  10,000  to  32,400.  Two  different  droplet  diameters  were  used  viz.,  60  fxm  and 
90  nm.  A  significant  range  in  the  Kolmogorov,  tm-bulent,  and  acceleration  Stokes  num¬ 
bers  was  covered.  The  times-of-flight  were  used  to  analyze  the  dispersion  measurements 
in  terms  of  true  Lagrangian  statistics.  It  was  fovmd  that  the  use  of  a  mean  time  of  flight 
to  present  the  data  incurred  quite  small  errors  relative  to  the  exact  Lagrangian  statistics. 
Dispersion  data  for  long  times  flight  over  the  range  of  nozzle  diameters,  particle  diame¬ 
ters  and  exit  velocities  were  analyzed  to  obtain  the  Lagrangian  particle  diffusivity.  The 
non-dimensional  diffiisivities  or  Peclet  niunbers  were  foimd  to  approeich  a  value  that  was 
similar  to  the  Eulerian  Peclet  nvunber  for  scalar  transport  in  the  jet.  Furthermore,  the 
dispersion  increased  hnearly  with  time.  The  results  suggested  that  some  economies  could 
be  achieved  in  the  simulation  of  dilute,  two  phase  turbulent  jets  and  similar  flows  over 
times  long  in  comparison  to  the  Lagrangian  integral  length  scale  by  taking  advantage  of 
their  self  preserving  nature. 


An  additional  goal  of  the  experiments  was  to  measure  the  dispersion  of  tagged  droplets 
in  a  Lagrangian  sense  within  a  turbulent  spray.  A  modified  experimental  apparatus  was 
designed  to  this  end.  IndiAddual  droplets  were  tagged  with  a  fluorescent  dye.  The  tagged 
droplets  were  injected  into  a  rotmd  turbulent  jet  into  which  a  dispersed,  liquid  spray 
had  been  added.  The  location  of  the  tagged  particle  was  followed  through  the  jet  with  an 
optical  detection  system  that  was  a  modification  of  an  earher  technique.  This  experimental 
approach  provided  a  simple  two  phase  flow  with  well  defined  boundary  conditions  without 
the  significant  compHcations  that  are  introduced  by  spray  atomization. 

Sprays  with  mass  loadings  from  0  up  to  7%  were  added  to  the  jet.  The  signal  to 
noise  ratio  under  these  conditions  was  about  110.  The  mean  square  displacement  of  the 
fluorescent  droplet  was  measured  at  varying  x/D  from  the  exit  of  the  nozzle.  It  was 
apparent  that  the  spray  contributes  to  the  breakdown  of  the  jet  with  larger  spray  loadings, 
resulting  in  greater  dispersion  of  the  tagged  fluorescent  droplets  that  were  added  to  the 
flow.  The  results  were  compared  with  the  predictions  of  a  Large  Eddy  Simulation  (LES) 
and  were  found  to  be  encomaging. 


□ 
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Theory  and  Computation. 

The  flow  fleld  in  turbulent  roimd  jets  and  the  paths  of  droplets  ranging  from  fluid  material 
points  to  heavy  particles  (hexEidecane  droplets)  were  simulated  numerically.  The  large  scale 
structmes  of  the  flow  fleld  were  resolved  and  the  scales  smaller  than  the  grid  were  modelled 
using  the  LES  approach.  The  particles  were  tracked  individually  in  the  Lagrangian  frame 
and  the  influence  of  the  particles  on  the  fluid  phase  was  represented  using  the  PIC  (particle 
source  in  cell)  method.  ,  ~ 


The  simulation  of  particle  laden  turbulent  jets  aUowed  the  computation  of  correlation 
Actions  for  pmticles  and  flmd  phase,  not  accessible  to  single  point  turbulence  closures. 
The  present  LES  model  for  the  turhulent  jet  produced  results  for  the  auto-correlatlon 
function  of  the  particle  velocities  in  good  agreement  with  measurements  for  the  case  of 
axial  velocity  components  shifted  in  time. 

The  numerical  simulation  of  dilute  sprays  was  performed  to  investigate  the  effect  of 
the  spray  on  the  near  field  of  a  turbulent  jet  at  a  moderate  Reynolds  number  Re  =  10  000. 
The  dispersion  of  tagged  particles  was  measured  and  simulated  numerically.  The  effect  of 
he  spray  on  the  flow  structure  emerged  as  reduction  in  the  size  of  the  largest  vortices  and 
an  increase  m  their  number  in  the  initial  region  of  the  jet. 

-nie  finite-difiFerence  method  for  the  simulation  of  the  fluid  phase  dynamics  was  im¬ 
proved  m  several  aspects.  The  original  Adams-Bashforth  time  integration  scheme  was 
second  order  accurate  and  unstable  in  the  inviscid  case.  Runge-Kutta  schemes  of  third 
and  fourth  order  accuracy  were  explicit  and  conditionally  stable.  Development  of  both 
v^^ts  was  begun  and  they  were  first  implemented  into  the  two-dimensional  Navier- 
Stokes  solvers  for  plane  and  axisymmetric  flows  with  excellent  results.  The  discretization 
of  the  convective  teiros  was  extended  to  windward  biased  schemes  up  to  seventh  order 
and  successfully  implemented  into  the  Runge-Kutta  versions  of  the  two-dimensional  and 
axisjunmetric  Navier-Stokes  solvers. 
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1.1  Computational  Approach. 

The  Navier-Stokes  equations  were  solved  using  an  accurate  finite-difference  method  dis¬ 
cussed  m  detail  in  Lienau  and  KoUmann  (1993)  and  the  final  report  (1989/92)  for  AFOSR 
rant  no.  89-0392.  The  void  fraction  for  dilute  sprays  was  approximated  by  unity  and  the 
equation  have  then  the  usual  form  for  incompressible  flows.  The  original 
method  proceeded  in  two  steps:  First  the  momentum  balances  were  solved  in  the  predictor 
step.  The  predicted  velocities  did  not  satisfy  mass  balance  and  a  second  step  was  necessary 
to  correct  velocity  ^d  pressure.  In  this  second  step  a  Poisson  equation  for  the  pressure 
correction  was  solved  and  velocity  and  pressure  were  updated. 

Momentim  and  mass  balances  were  set  up  in  cylindrical  coordinates  (r,  6,  C)  and  were 
discretized  using  ^curate  finite  difference  expressions.  Simple  grid  stretching  transfor¬ 
mations  were  apphed  to  the  r  =  r{f,/3)  and  C  =  C(C,/?)  directions  (/?  >  1  denotes  the 
stretching  parameter,  Lienau  and  KoUmann,  1993)  to  concentrate  the  grid  points  in  the 
region  of  interest  and  to  move  the  outer  boundary  as  far  out  as  possible  without  wasting 
too  many  ^d  points  The  time  derivative  was  discretized  using  the  second  order  accurate 
three  point  backward  (Adams  Bashforth)  scheme.  The  convective  terms  were  discretized 

mli  Tnd  M  ^  windward  biased  compact  differencing  scheme 

•  pressure  gradient  was  discretized  to  fourth  order  accuracy 

using  the  three  point  central  compact  scheme  (Lele,  1990)  and  the  viscous  terms  were 
discretized  to  fourth  order  accuracy  using  five  point  central  differencing. 

Mass  bailee  was  enforced  using  the  pressure  correction  method.  A  hybrid  method 
was  developed  for  the  solution  of  the  Poisson  equation  for  the  pressure  correction  that 
exploits  the  periodicity  of  the  flow  field  in  circumferential  ^-direction.  It  consisted  of  a 
our  h  order  compact  central  finite  difference  scheme  for  the  nonperiodic  dependence  on 
axial  and  ra^al  coordinates,  that  aUowed  the  appUcation  of  a  direct  solver,  and  a  Fourier 
spectral  method  m  the  periodic  0  direction.  The  direct  solver  had  the  advantage  that  it, 
provid^  the  solution  m  a  fixed  number  of  computational  steps  per  time  step  and  it  was 
more  efficient  for  large  problems  than  iterative  methods. 

Q  K  method  was  not  able  to  resolve  aU  scales  of  turbulence  and  a 

bub-Gnd-Scale  (SGS)  viscosity  model  was  added  to  represent  the  effect  of  the  unresol v(><i 
scales  of  motion  on  the  resolved  ones  (Reynolds,  1990).  The  Navier-Stokes  equations 


5 


were  written  for  filtered  variables  using  the  top  hat  filter.  The  closure  for  the  additional 
correlations  was  completed  by  neglecting  the  sum  of  Leonard  and  cross  stresses  and  the 
Smagorinsky  model  was  applied  to  the  SGS  stresses.  The  top  hat  filter  corresponded 
precisely  to  the  integration  for  the  force  term  representing  the  effect  of  the  particles  on  the 
fluid  phase  (see  section  1.1.2).  The  details  of  the  LES  approach  were  described  in  section 
1.1.1. 

Boundary  and  initial  conditions  were  critical  for  validity  of  the  simulation.  The  jet 
issuing  into  stagnant  surroundings  operated  in  an  inflnite  domain  with  the  boundary  condi¬ 
tion  that  velocity  was  integrable.  This  situation  was  approximated  by  a  cylindrical  domain 
with  radius  much  larger  than  the  jet  pipe  and  by  conditions  at  the  outer  boundary  that 
did  not  affect  the  flow  in  the  jet  significantly.  After  many  trials  a  condition  was  found 
that  worked  reasonably  well:  Set  the  normal  derivative  of  velocity  to  zero  on  the  outer 
part  of  the  boundary.  The  centerline  was  not  a  boundary,  since  a  point  on  r  =  0  was  inner 
point  of  the  flow  field.  The  finite  difference  grid  could  have  been  set  up  such  that  the  axis 
r  =  0  was  gridline,  which  required  a  detailed  analysis  of  the  equations  at  this  gridline  to 
remove  singularities  introduced  by  the  metric  of  the  coordinate  system,  or  such  that  r  =  0 
was  between  grid  surfaces.  The  latter  method  was  chosen  since  it  was  straightforward  to 
implement. 

The  pressure  gradient  in  the  direction  normal  to  the  computational  boundary  was 
set  to  zero  at  the  jet  inlet  and  exit  planes.  On  the  outer  boundary  the  pressure  and  the 
pressure  correction  were  set  to  zero.  This  provided  stability  for  the  jet  without  affecting  the 
solution  appreciably.  The  inlet  condition  for  the  jet  was  obtained  using  the  experimental 
results  of  Call  and  Kennedy  (1991)  or  Liepmann  (1991)  and  Liepmann  and  Gharib  (1992) 
and  the  entire  flow  field  was  initialized  with  the  inlet  profile,  representing  an  inviscid  flow 
initially.  This  type  of  initizd  condition  weis  superior  to  starting  the  jet  flow  by  increasing 
the  mass  flow  rate  imtil  the  desired  steady  state  value  was  reached.  The  flow  set  up 
this  way  was  first  dominated  by  the  start  up  vortex,  which  moves  through  the  flow  field 
creating  an  axial  flow  similar  to  the  inviscid  initial  condition,  but  requires  significantly 
more  CPU-time  without  significant  change  of  the  later  flow  development. 

Disturbances  were  introduced  at  the  inlet  to  initiate  and  accelerate  the  roll-up  process 
of  the  cylindrical  shear  layer.  The  distm-bance  V(^{0,t)  was  chosen  such  that  it  was  only 
a  fimction  of  9  and  time  similar  to  Rai  and  Moin  (1991)  and  satisfies,  therefore,  mass 
balance.  It  was  given  by 


Ai^nsin{k9  +  <l>e)sin{ — —  +  (})t) 

k=l  n+l 


(1) 


where  =  cvjetTjet^^  ^  =  0-25,  7  =  0.5,  c  =  0.1  Two  modes  in  both  9  and 

time  were  selected  (i.e.  K  =  2,  N  =  2). 

The  analysis  for  the  Adams  Bashforth  time  integration  method  and  the  compact  four 
point  windward  biased  scheme  for  the  convective  terms  revealed  that  the  maximum  stable 
Courant  number  was  0.5.  Combining  the  time  step  restrictions  for  convective  and  viscous 
terms  lead  to  the  maximal  time  step 


At  <  0.5 


min 

X 


dr  \Vr 
dr  Af 


1  H 

rA9 


^Cl^d  I  1  (fdr^  1 
dC  AC  3Rey  dr  ^  Af2 


which  was  used  to  control  the  time  step. 


1 

r^A9‘^ 


(2) 
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1.1.1  Large  Eddy  Simulation 

The  Large  Eddy  Simulation  of  the  turbulent  flow  in  round  jets  was  the  most  realistic 
approach  for  the  prediction  of  the  statistical  properties  of  the  flow  field  and  the  simulation 
of  Lagrangian  particle  dynamics  in  this  type  of  flow,  since  direct  simulation  of  turbulence 
in  roimd  jets  at  high  Reynolds  numbers  was  still  not  within  the  capabilities  of  present 
computers.  The  simulation  of  turbulent  flow  fields  in  roimd  jets  in  the  present  project  was 
based  on  highly  accurate  finite-difference  methods,  which  offer  the  flexibility  necessary  for 
the  treatment  of  non-periodic  jet  flows  emitting  from  nozzles  and  the  consideration  of  a 
variety  of  exit  geometries  and  velocity  profiles.  Furthermore,  the  influence  of  disturbances 
created  at  the  jet  pipe  exit  on  the  flow  development  was  able  to  be  studied  in  detail.  This 
aspect  of  the  project  was  described  in  detail  in  the  final  report  for  grant  AFOSR  89-0392 
and  further  details  can  be  found  in  Lienau,  Kennedy  and  Kollmann  (1993)  and  Lienau  and 
KoUmann  (1993)  of  the  publication  list. 

The  accomplishments  of  the  current  review  period  covered  the  implementation  of  the 
LES  model  and  improvements  of  the  solution  procedure.  The  numerical  treatment  of  the 
coordinate  axis  r  =  0  for  unsteady  flows  without  symmetries  was  analyzed  (details  in  the 
appendix  of  the  annu8il  report  for  1992/3)  and  a  satisfactory  method  was  found  to  avoid  the 
loss  of  accuracy  near  the  axis.  The  second  and  main  contribution  was  the  implementation 
of  a  LES-model  for  the  non-resolved  scales  of  the  turbulence.  The  Navier-Stokes  equations 
axe  written  for  filtered  variables 


dx'G(x  —  x',t) /(x',  t) 


(3) 


where  G{x  —  x',  t)  denotes  the  filter  function  and  /(x,  t)  a  dependent  variable.  The  Navier- 
Stokes  equations  (in  Centesian  coordinates  for  convenience)  were  filtered  emd  emerge  as 


dvg 

dXg 


=  0 


(4) 


and 


,  d  r-  .  1  ap  .  d^Vg  d 

dt  dxp  “  ^  pdxg  ^  dx^dxp  5x^ 


{Lg^  +  Cot^  -1-  Rg^) 


(5) 


The  modified  filtered  pressure  was  defined  by 


P=P+  oP^apVgV'ff 


(6) 


and  the  correlations  of  filtered  and  sub-grid-scale  motion  (v'^  =  Vg  —  Vg)  were  given  by 


■^a/3  — 

(7) 

Cgfi  =  V'J}0  -t-  VgV'p 

(8) 

-a/3  =  -  ^8a0V^V>^ 

(9) 

The  present  filter  was  the  top  hat  filter  and  the  closure  for  the  correlations  is  given  by 

Pa/3  +  Cg^=0  (10) 
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and  the  Smagorinsky  model  for  (Reynolds,  1990) 

Ra/3—  ~  2i/7’S'ay(3 

where  the  eddy-viscosity  was  defined  by 


and 


and 


=  (cAf(2S^0S^^)i 
A  =  (AxAyAz)^ 


(12) 


(13) 


>a/3 


2\dx^  dxa) 


(14) 

was  the  rate  of  strain  produced  by  the  filtered  velocity  field 
the  Iw  o^hf  Md"'nhai'^a  i"’’  »  ‘he  structure  of 

b“d°  Ket:“  ;:T  £r  « 

•  sL  TT  y  V  yyij  at  the  nominal  Reynolds  number  of  i?p  =  !'=»  nnn  ^ 

£s„ 

t“ft;m  Jr  02Tf  -or-plfyslhen 

a  noticeable  role  for  dispeSr^  5  muL^at^xTef 

the  experimental  (mean)  time  of  flight  was  mnrh  1  f  fk  ^  spersion  calculated  with 
of  th#>  T  T?Q  vrJe  •\  {i.  much  closer  to  the  measurements.  The  effect 

the  LES-viscosity  on  the  dispersion  was  shown  in  Figure  6  Thp  in^vroco/a  • 
due  to  the  LES-model  (101  decrej^Pd  thA  d?=rv  •  J-he  increase  in  viscosity 

measurements.  Fi^X!  it  wt  fou^^^^^^  fheTan 

fairly  independent  of  the  method  of  calculating  the  ttoe  of  %ht!^  “ 
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1.1.2  Spray  simulations. 

The  turbulent  jet  used  in  the  experiment  for  the  study  of  droplet  dispersion  had  a  suffi¬ 
ciently  low  Mach  number  such  that  the  flow  of  an  incompressible  fluid  containing  particles 
was  considered  incompressible  for  the  numerical  simulations.  The  Navier-Stokes  equations 
for  the  incompressible  fluid  phase  were  set  up  for  the  conditions  of  a  flow  laden  with  fine 
particles.  The  volvune  occupied  by  the  particles  per  unit  volume  denoted  by  0{x,t)  was 
defined  by 

=  (15) 

where  V  was  the  volume  of  a  subdomain  V  of  the  flow  field,  that  was  contracted  to  the 
limit  Vq  centered  at  the  point  x,  and  Vp  was  the  voliune  of  all  particles  Vi  in  V.  The 
limit  value  Vq  for  the  subdomain  V  in  (5)  for  dilute  sprays  of  particles  much  smaller  than 
the  smallest  length  scale  of  the  flow  can  be  chosen  sufficiently  small,  such  that  the  void 
fraction  can  be  regarded  as  local  variable.  The  limit  was  then  well  defined  and  the  void 
fraction  emerged  for  spherical  particles  in  the  form  (Dukowicz,  1980) 


N{x) 

e{x,t)  =  i- 

i=l 


47r  r? 


(16) 


where  the  smn  extends  over  all  particles  in  Vq  and  N{x)  was  the  nxmiber  of  particles  in  Vq 
centered  at  x.  The  mass  balance  for  the  two-phase  system  had  then  the  form  (Dukowicz, 
1980) 


d 

dxa 


(H)  =  o 


and  momentum  balance  was  given  by 


dt  ^dXa 


1  dp 
pf  dXa 


+  i 


where  the  superscript  /  indicated  the  fluid  phase.  The  last  term  in  the  momentum  balance 
represented  the  effect  of  the  particles  on  the  fluid  momentiun  as  defined  below.  The  case  of 
dilute  sprays  of  small  particles  was  considered  in  the  following,  where  the  approximations 
were  made:  the  particle  size  was  much  smaller  than  the  smallest  length  scale  of  the  fiow 
and  that  9  =  1.0  holds.  Hence,  the  displacement  effect  of  the  particles  was  neglected,  but 
not  the  forces  exerted  by  the  particles  on  the  fluid.  The  Navier-Stokes  equations  emerged 
then  in  the  form  (Elghobashi  and  Truesdell,  1991,1993,  Eaton,  1994) 


d 

dXa 


(vl)  =  0 


(17) 


Momentum  baleince 


+  V' 


M 

^dx„ 


1  dp 
pf  dXa 


dxp^"'^ 


which  contained  the  effect  Ma  of  the  particles  on  the  fluid  phase. 


(18) 
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Particle  dynamics. 

Solid  spherical  particles  in  a  gaseous  fluid  phase  were  considered.  The  position  of  particle 
p  at  time  t  was  denoted  with  XFify  and  its  velocity  with  The  position  and  velocity 

of  a  particle  were  related  by  definition 

dXP 

-^(t)  =  (19) 

with  the  initial  condition  2C^(0)  =  (x  denoted  the  release  position  of  particle  p).  The 
particle  velocity  vP(t)  was  different  from  the  Eulerian  velocity  field  v/{x,t)  of  the  fluid 
phase,  which  was  the  solution  of  the  Navier-Stokes  equations  (17),  (18)  obtained  with  the 
finite-difference  method  described  in  a  later  section.  The  dynamic  equation  determining 
the  location  of  particles  with  nonzero  mass  was  a  consequence  of  Newton’s  second  law  with 
approximations  for  the  forces  acting  on  the  particle  (see  Odar  and  Hamilton,  1964).  The 
equation  for  the  peirticle  velocity  contained  terms  representing  drag  and  virtual  mass,  the 
Basset  term  and  buoyancy. 


2  pp  Dtf 


pg  DyJ 
Pp  Dif 


dr  D 


(t-r) 


^  DXr 


(„/_„P)  +  (££_l)g 

Pp 


(20) 


The  subscript  /  denoted  fluid  (gas)  and  p  particle  properties,  the  substantial  derivatives 
were  to  be  taken  accordingly  as  following  a  material  point  of  the  fluid  or  the  particle,  d 
was  the  particle  diameter,  u  was  the  gas  viscosity  and  uf,  pf  were  taken  as  the  values 
undisturbed  by  the  particle  at  the  location  of  the  particle.  The  coefladents  Cd,  Cj  and 
Cb  were  in  genered  functions  of  the  particle  Reynolds  number 


Rep 


(21) 


and  the  acceleration  number  (see  Faeth,  1983  and  Hansell  et  ed.,  1992) 


^\uf  -uP\ 


\y/  —  uPp 

The  drag  coefficient  was  given  by  (Putnam,  1961) 


24  Ttp^ 

Cd  =  +  for  Hep  <1000 


(22) 


(23) 


The  present  calculation  of  the  particle  dynamics  took  only  drag,  virtual  mass  and  buoyancy 
into  account.  The  Basset  term  was  rather  expensive  to  calculate  and  will  be  implement  t'<l 
at  a  later  stage.  Therefore,  particle  velocity  was  the  solution  of 


dup 

dt 


3  Cp  Pg 

^  d  Pp 


1Lg-Up\{Ug-Up)^- 


Pp  Dtg 


(21) 
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For  the  present  computations  the  properties  for  rf  =  36  /xM  hquid  hexadecane  particles 
were  used,  which  were  given  by  pp  =  775^.  The  relevant  properties  of  air  were  u  = 
1.45(10)-=!^,  ^,  =  1.19^. 

The  time  integration  for  (19)  and  (20)  was  carried  out  simultaneously  with  the  inte¬ 
gration  of  the  Navier-Stokes  equations.  The  second  order  accurate  improved  Euler  method 
was  applied  to  the  time  integration  of  the  particle  equations,  which  was  the  same  order 
of  accuracy  as  the  current  time  integration  scheme  for  the  Navier-Stokes  equations.  The 
values  of  velocity  at  an  arbitrary  location  inside  a  grid  cell  were  required  for  this  integra¬ 
tion  method.  They  were  interpolated  using  a  second  order  polynomial  (Bezier  Spline)  in 
three  dimensions.  This  involved  using  a  three  point  stencil  consisting  of  the  two  points 
defining  the  cell  containing  the  particle  location  and  the  point  ahead  or  behind  this  cell 
depending  which  point  was  closest.  The  three  dimensional  interpolation  scheme  was  rather 
complicated  and  only  the  overall  Bezier  spline  in  one  dimension  will  be  given 

=  7  +  2a6V'j  +  b\+i)  (25) 

where  'll;  was  the  variable  being  interpolated,  a  was  the  distance  from  the  particle  to  point 
j  +  1  and  b  was  the  distance  from  the  particle  to  the  point  j  —  1  and  c  was  the  distance 
from  point  j  —  1  to  j  -f- 1  and  the  particle  is  between  j  —  1  and  j. 


Interaction  of  fluid  and  particles 

The  particle  dynamics  were  driven  by  the  forces  exerted  by  the  fluid  phase  on  the  particles. 
The  effect  of  the  particles  on  the  fluid  phase  was  represented  by  the  force  term  M  in  the 
momentum  balance  (18)  and  the  displacement  effect  of  the  particles  described  by  the  void 
fraction.  If  the  momentum  balance  (24)  for  the  particles  was  written  as 


dt 


=  FP 


(26) 


where  was  the  siim  of  all  forces  acting  on  a  particle,  then  it  can  be  deduced  from  the 
spray  equation  (see  Williams,  ch.ll,  1985),  that  the  fluid  phase  was  exposed  to  the  force 


Ma  =  I J  dvr^F^  f{r,x,v,t) 


(27) 


exerted  by  the  particles.  Here  f{r,x,v,t)  was  the  probable  number  of  particles  in  the 
radius  range  dr  about  r  located  in  the  spatial  dxa  range  about  Xa  with  particle  velocities 
in  the  range  dva  about  Va  satisfying 


N  = 


dv  f{r,x,v,  t) 


where  N  was  the  total  number  of  particles  in  the  domain  of  integration  V  with  volume 
y  >  0.  The  solution  f{r,x,v,t)  of  the  spray  equation  was  in  general  not  known  and 
numerical  simulation  methods  were  applied  to  obtain  the  information  on  the  force  term 
(27).  In  a  direct  simulation  method  all  particles  would  have  been  tracked  individually  and 
the  distribution  function  / (r,  x,  v,  t)  would  be  approximated  by 


f{'f'-,x,v,t)  =  y^^<^(r-r*)S(x-x*)S(v  —  v*) 

i=l 
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force  A/,  was  integrated  over  it.  It  foUowed  that  selected  and  the 


^  fdxM^  =  -^f'YLF< 
PlkVc'‘ 


(28) 

holds,  where  the  £<  were  the  forces  on  the  K  individual  particles  in  the  grid  ceil  V  with 
ume  ir  The  particle  volumes  were  for  spherical  particles  given  hy  K  =  This 

w:rc:,i  rthr  “>■  p-S- 

nnf  numerical  simulations  at  the  Reynolds  number  (10,000)  of  the  experiments  could 
not  have  resolved  the  complete  range  of  scales  present  in  the  turbulent  flow  aTthe  LES 
model  descnbed  in  section  1.1.1  was  appUed  to  simulate  the  effect  of  the  unresolved  part  of 

T  PQ  •  J  However,  at  this  level  of  mass  loading  the  dispersion  produced  bv  the 
LES  simulation  m  Figure  7  showed  clearly  the  phenomenon  observed  in  the  experiiints 

J/D  -  nwt’’  Tt.'"  further  downstream  (at  about 

it 1  ’I-  “v  «*P«"nient,  but  the  shape  was  in  good  agreement  with 
.  The  explanation  for  this  phenomenon  can  be  found  if  the  detailed  flow  structure  in 
his  region  was  considered.  The  particles  were  fed  into  thp  finw  of  fVi  + 
uniformly.  Little  change  of  the  uniformity  of  the  parti^k  rcToTs  "ctn 

of  the  jet  was  observed  for  the  first  five  diameters.  Once  the  roll  im^f  M.  ! 

th?  interact  the  particles  were  transported  SaUy  ^d 

the  particle  number  density  becomes  increasingly  non-uniform.  At^the  location  of  the 
dispersion  plateau  seen  in  Figure  7  in  the  ranee  10  <  t/D  <rifi  •  «  the 

the  particle  di^itrihnfinr,  rvff  il  •  ^  13  a  sigmficant  thinning  of 

the  Kalvin-Hcli^oltz  instabiUty  of  the  cyhndricaTshir  iTyer  a  m^ti  Ssfswhtf  the 

jet  a^s.  Inspection  of  the  velocity  field  showed  that  this  was  indeed  the  case  The 
flow  accelerated  and  moved  towards  the  jet  axis  as  indicated  in  the  particle  distribution 
f  several  time  instances  were  analyzed  it  became  evident  that  the  statistical  result  for 
he  p^icle  dispersion  was  the  plateau  in  the  dispersion  emerging  in  Figure  7  Figure  8 
showed  two  scatter  plots  of  the  numerically  simulated  distribution  of  thf  particles  It  the 
^al  statioi^  x/D  =  20  and  x/D  =  28  downstream  of  the  plateau  region  "rel  in 

^ompoueut  of  voctici^,  which  i,  responoiblf^tl:;!;"  ut^Tw 

®,7u  V“  j®*-  ^'Sure  9  shows  the  case  with  partiLs 

Su?  p  Jri^  *^®  conditions 

hout  particle  injection.  It  was  seen  that  the  presence  of  particles  led  to  more  ring 

plftidl  "®^®"  “  the  flow  withou? 

The  flow  structures  in  the  full  computational  region  1  <  x/D  <30  were  presented 
as  iso-hnes  of  the  axial  and  radial  velocity  components  in  Figure  11  without  bi^kscatter 
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effect  (the  term  (28)  was  not  included  in  the  momentmn  balance)  and  Figure  12  with  this 
effect.  The  reduction  in  scale  was  clearly  visible  but  the  backscatter  effect  was  noticeably 
smaller  in  the  downstream  region.  The  presence  of  the  particles  induced  a  larger  number 
of  ring  vortices,  which  were  smaller  and  thus  spaced  more  closely  than  in  the  flow  without 
the  backscatter  effect. 

1.1.3  Time  integration  schemes. 

The  second  order  accurate  explicit  Adams-Bashforth  time  integration  scheme  was  known 
to  be  unstable  for  inviscid  flows  and  thus  not  suitable  for  long  time  integration.  It  was, 
therefore,  necessary  to  replace  this  scheme  with  schemes  that  were  at  least  conditionally 
stable  and  at  least  as  accurate.  Several  single  and  multi-step  time  integration  schemes 
were  analyzed  (details  can  be  foimd  in  the  appendix)  and  two  Runge-Kutta-type  schemes 
were  successfully  tested.  The  fourth  order  Runge-Kutta  scheme  was  described  in  detail  in 
the  appendix. 

A  standard  third  order  Runge-Kutta  method  for  an  ode  y'  =  /(x,  y)  could  be  set  up 
as  follows 

=  y”  -f  7oA:o  +  ')\ki  -t-  72^2 
where  the  increments  were  defined  by 


ka  =  Atf{y'^) 

ki  =  At/(j/"  +  Q!ifco) 

^2  =  Atf{y'^  +  f32ikQ  +/?22fci) 
and  the  constraints  for  the  coefficients 


/?22  =Ot2—  /?2l 


70  +  7i  +  72  =  1 
7i«i  +720:2  =  ^ 

7iO:i  +  720:2  =  ^ 

72P22OC1  =  i 

There  was  no  unique  set  of  constants  and  several  variants  of  the  third  order  method  were 
available.  The  application  of  a  Runge-Kutta  time  integration  scheme  to  the  solution  of  the 
Navier-Stokes  equations  in  three  dimensions  introduced  a  new  aspect  to  the  construction  of 
such  a  scheme:  Minimization  of  storage  requirements.  A  possible  solution  to  the  constraint 
equation  was  then  given  by 


—  2’  ^21  =  2’  /^22  =  — - 

2  5 

70  =  3,  71  =  3. 
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4 

72=-3 


The  procedure  was  tested  for  the  solution  of  Biurgers  nonlinear  second  order  pde  with  sixth 
order  compact  finite  difference  approximations  for  the  convective  and  diffusive  terms  and 
the  present  third  order  Runge-Kutta  time  integration  scheme.  Absolute  error  bounds  were 
established  since  the  analytic  solution  for  the  periodic  initial  value  problem  was  known 


u{x,t) 


2  M(x,t) 
Re  N(x,t) 


where 


M(x,t)  =  sin(x)exp(-^)  +  isin(2x)exp(-^)  +  ^sin(3a:)exp(-^) 

+-  sm(4x)  exp(-— )  +  -  sin(5x)  exp(-— ) 

Ar(x,t)  =  2  +  cos(x)exp(-4-)  +  7Cos(2x)exp(-4^)  +  ^  cos(3x)  exp(--^) 

xie  4  He  2  Re 

cos(4x)  exp(-^)  +  i  cos(5x)  exp(-^) 

Restricting  the  time  step  by 

At  <  Re^ 

the  error  was  computed.  Selecting  Re  =  100  and  the  time  t  =  20  the  results  were  shown  in 
Figure  13  and  Figure  14  proving  that  the  time  integration  weis  indeed  third  order  and  the 
spatial  discretization  was  sixth  order  as  claimed.  The  fourth  order  Rimge-Kutta  method 
for  the  Navier-Stokes  equations  was  discussed  in  detail  in  the  appendix. 

The  application  of  the  Runge-Kutta  methods  to  the  Navier-Stokes  equations  was  first 
carried  out  first  for  the  case  of  two-dimensional  flows  in  Cartesian  and  axi-symmetric 
coordinates  in  order  to  save  CPU-time  and  speed  up  the  testing.  The  case  of  third  order 
Runge-Kutta  time  integration  and  compact  spatial  discretization  to  fourth  order  for  the 
viscous  terms  and  the  pressure  gradient  and  seventh  order  windward  biased  scheme  for  the 
convective  terms  applied  to  the  axis-symmetric  round  jet  was  presented  in  Figure  15.  The 
Reynolds  number  was  Re  =  5000  and  the  time  t  =  46  dimensionless  units.  The  absolute 
value  of  vorticity  was  plotted  for  the  first  fifteen  diameters  of  the  jet.  No  disturbances 
were  introduced  and  the  natural  instability  of  the  axi-symmetric  shear  layer  produces 
vortex  roll-up  and  vortex  merging.  The  maximum  norm  of  vorticity  in  Figure  16  shows 
the  variation  with  time  due  to  the  roll-up  and  merging  processes.  The  check  on  mass 
balance  was  presented  in  Figme  17  in  terms  of  the  integral  of  the  divergence  of  velocity  as 
broken  line  (the  exact  value  is  zero)  and  the  L2  of  the  divergence.  It  was  seen  that  both 
are  satisfactory.  The  fourth  order  Rimge-Kutta  time  integration  method  was  implemented 
into  the  Ceirtesian  version  of  the  Navier-Stokes  solver.  Some  results  for  two-dimensional 
flows  were  shown  in  Figure  18  and  Figure  19.  The  vector  plot  for  velocity  for  the  plane 
jet  in  Figure  18  at  the  Reynolds  number  Re  =  2500  and  t  =  25  indicated  the  difference 
between  axi-symmetric  and  plane  jets  as  the  latter  has  a  strong  tendency  to  meander.  The 
last  test  example  was  the  backward  facing  step,  which  was  simulated  for  the  Reynolds 
number  (based  on  the  step  height)  Re  =  5000  and  the  ratio  of  step  height  to  channel 
width  of  one  half.  The  flow  was  simvdated  starting  with  zero  velocity  and  increasing  the 
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inflow  velocity  to  its  steady  state  value  in  two  time  units.  Two  times  were  shown  in  Figure 
19,  t  =  15  in  the  upper  graph  and  t  =  25  in  the  lower  one.  The  flow  was  not  yet  fully 
developed  as  oscillating  separating  vortices  dominate  the  flow.  It  was  concluded  from  the 
results  presented,  that  the  Runge-Kutta  time  integration  schemes  were  very  promising  for 
the  full  three-dimensional  simulations.  These  schemes  were  implemented  into  the  fully 
three-dimensional  solver. 

1.2  Experimental  Approach. 

The  experimental  method  has  built  upon  the  pioneering  work  of  Call  and  Kennedy  (1991) 
wherein  a  method  was  developed  and  reported  for  the  Lagrangian  tracking  of  particles  in 
a  turbulent  jet.  Turbulent  jets  pose  a  particular  problem  for  the  description  of  particle 
dispersion  in  that  the  mean  flow  field  undergoes  a  rapid  change.  Hence,  unresponsive 
particles  may  find  themselves  lagging  the  development  of  the  flow  field.  Although  there 
have  been  a  number  of  calculations  of  the  dispersion  of  particles  in  turbulent  jets  with 
Lagrangian  tracking  methods,  there  have  been  few  experimental  studies  of  a  Lagrangian 
nature.  Call  and  Kennedy  (1994)  used  the  same  experimental  apparatus  as  the  one  that 
was  used  during  this  reporting  period  to  investigate  Lagrangian  dispersion  in  a  turbulent 
shear  flow.  Various  sizes  of  vaporizing  and  non-vaporizing  droplets  were  used.  However, 
the  experiments  employed  a  limited  range  of  particle  and  fluid  time  scales.  The  new 
experiments  covered  a  wider  range  of  time  scales  by  variations  in  the  particle  size,  the 
nozzle  diameter  and  the  jet  exit  velocity. 

The  experiment  used  hexadecane  droplets  that  were  more  than  1000  droplet  diameters 
apart.  Consequently,  there  were  no  pgirticle  interactions.  The  peirticles  were  non  vaporizing 
because  the  air  was  at  room  temperatme  (23  ®C)  and  the  particles  were  measured  in  the 
flow  for  a  relatively  short  time  (<^  1  s).  Monodisperse  droplets  were  produced  by  a  piezo¬ 
electric  crystal.  The  crystal  was  a  hollow  cylinder  with  a  fuel  inlet  on  one  side  and  a  nozzle 
on  the  other.  A  pulse  generator  supplied  a  pulse  of  variable  voltage,  frequency  and  width 
to  the  crystal.  This  voltage  caused  the  crystal  to  contract  and  push  a  small  amount  of 
liquid  through  the  nozzle.  After  the  droplet  was  created,  it  was  allowed  to  fall  freely  in 
the  droplet  shroud.  The  shroud  ensured  that  the  droplet  remained  on  the  centerline  of 
the  flow.  The  air  was  straightened  through  screens  and  a  honeycomb  before  it  came  into 
contact  with  the  droplet.  The  air  and  droplet  were  then  accelerated  through  the  converging 
nozzle  and  the  droplet  entered  the  jet  on  the  centerline. 

The  droplet  detection  system  used  a  He-Ne  laser  that  crossed  the  centerline  of  the 
nozzle  exit.  As  a  droplet  left  the  nozzle,  it  crossed  the  laser  beam,  scattering  light  that 
was  focused  onto  a  photo  diode.  This  voltage  signal  was  used  to  trigger  the  data  collection 
system.  Two  laser  sheets  were  formed  by  em  zirgon-ion  laser  with  associated  optics.  As  the 
droplet  passed  through  the  sheets,  a  position  sensitive  photo  multiplier  tube  detected  the 
scattered  light.  The  photomultiplier  tube  had  four  anode  outputs  whose  magnitudes  were 
proportional  to  the  proximity  of  the  centroid  of  the  scattering  image  to  the  corresponding 
side  of  the  detector.  The  signals  were  amplified  before  being  digitized.  At  least  10(K) 
droplets  were  measured  in  order  to  give  statistically  meaningful  results.  The  position  of  a 
particle  was  determined  from  scattering  from  one  laser  sheet. 

Measurement  of  the  axial  component  of  the  particle  velocity,  tt,  required  the  use  of 
two  laser  sheets.  The  distance  between  the  center  of  the  sheets  was  measured  to  give  tli<‘ 
length  AL„.  The  time  between  the  peaks  of  scattering  signals  was  calculated  as  follows. 
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Table  1.  Particle  Peclet  Numbers  for  60  fim  Droplets 
D  =  7  ixm 

Rej  =  20, 000  Rsj  =  30, 000 

39  43 


D  =  12.6  /im 

Rej  =  20, 000 

53 


Rej  =  30, 000 
59 


A  Gaussian  curve  fit  was  applied  to  each  signal  with  a  linear  regression  performed  on  31 
data  points.  This  process  calculated  the  time  of  flight,  At^,  between  the  sheets  to  A  1/4 
of  one  clock  pulse.  Hence,  the  axial  velocity,  u  ,  was  found  from: 

u{xi)  =  ALu 

Ail 

In  principle,  the  system  should  be  capable  of  measuring  the  radied  component  of  the  droplet 
velocity.  This  component  of  velocity  is  more  relevant  to  the  radial  dispersion  than  the  axial 
velocity.  However,  it  was  found  that  in  practice  it  was  difiicult  to  resolve  accmately  the 
rather  small  velocities  in  this  direction  with  the  present  apparatus. 

A  TSI  model  1053B  hot  wire  anemometer  was  used  for  the  velocity  measurements  of 
the  air  jet.  The  exit  profiles  for  both  the  7  mm  cuid  12.6  mm  nozzles  exhibited  a  uniform 
exit  velocity  (top-hat)  from  both  nozzles.  The  radial  distance  X2  is  normalized  by  the 
nozzle  diameter  D  ;  the  velocity  u,  is  normalized  by  the  centerline  exit  velocity  Uc- 

The  axial  mean  velocities  and  turbulence  intensities  within  the  air  jet  were  also  mea¬ 
sured.  It  was  apparent  that  the  jet  at  a  Reynolds  number  of  10,000  may  not  have  been 
fully  developed  while  the  flow  at  the  other  Reynolds  numbers  appeared  to  exhibit  the 
characteristics  of  high  Reynolds  number,  self  preserving  turbulent  flow.  Measurements  of 
the  radial  mean  velocity  profiles  of  the  air  at  different  axial  positions  (not  shown)  indicated 
that  they  were  self-similar. 

Measurements  of  particle  statistics  were  obtained  for  a  range  of  characteristic  particle 
response  times,  nozzle  diameters  and  jet  exit  velocities.  It  was  foimd  that  increasing 
the  Reynolds  number  of  the  jet  by  increasing  the  exit  velocity  resulted  in  less  particle 
dispersion  in  the  Eulerian  frame  as  a  result  of  the  reduced  response  of  the  particles  to  the 
turbulence  and  the  reduced  time  for  a  particle  to  travel  from  the  nozzle  exit  to  a  given 
location  downstream.  Conversely,  the  particle  dispersion  in  the  Lagrangian  frame  i.e.,  as 
a  function  of  time  of  flight,  exhibited  an  increase  with  Reynolds  number.  Dispersion  as 
a  function  of  time  of  flight  was  quadratic  for  short  times  of  flight.  The  function  became 
linear  for  longer  times.  It  was  argued  that  when  the  local  particle  Stokes  number  was 
less  than  unity,  the  particle  would  behave  like  a  fluid  particle.  The  hnear  function  was 
an  indication  of  the  plausibility  of  the  assumption  of  Batchelor  (1957)  and  of  Monin  and 
Yaglom  (1979)  that  the  Lagrangian  statistics  in  a  free  shear  flow  would  be  self-preserving  in 
the  same  manner  as  the  Eulerian  statistics.  The  Lagrangian  Peclet  numbers  of  the  particles 
(shown  in  Table  1)  approached  the  estimated  v^ue  of  the  Eulerian  Peclet  number  for  a 
scalar  in  a  roimd  jet  i.  e.,  aroimd  50.  The  results  suggest  that  the  calculation  of  particle 
dispersion  in  self-preserving  flows  such  as  jets  can  be  simplified  considerably  once  the  local 
particle  Stokes  nxmiber  is  close  to  one.  The  simplification  could  be  achieved  by  adopting 
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a  constant  particle  difiusivity  or  equivalentlv  in  fho  i?  iu 

that  the  dispersion  of  a  given  size  class  nf  r»  +•  i  •  ^  ^  round  jet,  by  assuming 

the  dispersion  is  equivalent  to  the  variance  of ^  Gauss'*^^^^^  k  Then 

particle  location.  Hence,  the  distribution  of  partfcir^Toru^ fl  ’’  u’' 

descnbed  as  a  function  of  the  mean  timp  nf  fl-  v,+  I 

cane,  it  would  not  be  necessary  to  undertake  thi  rim"*’  location.  In  that 

Mood^^LlrSr"  within  the  paper  by 

1*2.1  Dispersion  in  a  Spray* 

c^cltd"^?;;  ‘»e  lifetoe  of  the  grant  was 

Fluorescent  tagging  was  Ised  pSk  ^  turbulent  spray. 

An  ultrasonic  atomizer  was  used  orodnnp^  from  the  surrounding  droplets. 

(see  Figure  21).  The  spray  was  introdnJlH  t  ^  5^u  around  50  ixm 

(Pigurf  22).  4e  chXrZS^oZTr^  ‘“n''*  ®7  “  -*»“ber 

through  the  central  tube  to  prevent  the  snrav  hitt^  ^  through  the  walls  and 

considerable  amount  of  work  was  necessary  in  order^f  creating  drips.  A 

the  spray  to  be  introduced  without  dripniL  frZ^Th  that  permitted 

seeded  with  fluorescein  dye  were  created^at^hp  t’  droplets  that  were 

droplet  generator.  The  fluLIce^Tolt  n  the  chamber  with  a  piezoelectric 

the  centerline.  As  it  passed  through  an  aJT  ^  through  the  spray  to  leave  the  nozzle  on 
a  longer  wavelength.  The  fluorescence  was  detecT^  ®”“*ted  fluorescence  at 

tube  to  yield  a  Measure  of  fh/SraTdt^^  photonrultiplier 

strong  background  Mkl^'teringTom  ’sS  Th”  ‘It®  ''cy 

with  a  Raman  filter  followed  by  f  colored  T  *^^^^ground  scattering  was  blocked 

offer^  up  to  6  orders  of  nn^^r S^ft^  C  ^::ieSr““ 

thesoaZrrrs^a^ilrnri^^^^^ 

droplet  velocity.  Higher  mj 

capability  to  form  denser  sprays.  ^  ^  proceeding  on  improving  our 

Sprays  with  mass  loadings  from  0  un  to  7%  tiawra  k  j  j  i  .  . 

noise  ratio  under  these  conditions  was  about  110  Tb  the  jet.  The  signal  to 

fluorescent  droplet  was  measured  at  varvine  v /n  f  fk  square  displacement  of  the 
is  shown  in  Fi^e  7  for  aCe  air  ilZn  f  ^  dispersion 

apparent  that  the  spray  contributes  to  the  bm'akdT^^orthe  loadings.  It  was 

resulting  in  greater  dispersion  of  thp  taa,y  ri  a  ^  ^  ^  loadings, 

flow.  Anoth»  interesting  Sru^  ‘f  ^dded  to  th 

x/D  «  8  -  10.  *"  ‘he  dispersion  at  around 
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Appendix:  Analysis  of  time  integration  schemes 

Single  step  and  multi-step  time  integration  schemes  were  considered.  Accmacy,  stability, 
number  of  operations  per  step  and  storage  requirements  were  the  relevant  properties  for 
the  application  to  the  solution  of  the  Navier-Stokes  equations.  The  general  theory  for  the 
numerical  solution  of  systems  of  odes  using  single  step  methods  were  found  in  Grigorieff 
(1977a)  and  using  multi-step  methods  in  Grigorieff  (1977b). 


Adams- Bash  forth  multi-stev  explicit  schemes. 
The  ordinary  differential  equation 


du  ,,  , 


with  the  initial  condition  w(0)  =  mq  was  considered.  Assmning  that  at  least  m  >  0 
steps  have  been  computed,  m  values  Uj^k  =  u{tj+k)  and  fj+k  =  f{tj+k,Uj+k)  for  k  = 
0,  •  •  • ,  m  - 1  were  used  to  construct  an  explicit  scheme  for  the  calculation  of  Uj+rn-  The  ode 
was  integrated  and  the  right  side  f{t,  u)  was  approximated  by  the  interpolation  polynomial 
P{t)  through  the  points  {tj,  /j),  •  •  • ,  {tj+m,Uj+m) 


'^j+tn  —  '^j+m—1  "b  J  dt'  P{t'') 

The  polynomial  P(t)  was  given  by  Newton’s  formula  using  backward  differences  as 

m-l 


where 

and  VVi  =  fi 
The  integration  led  to 


X  X 


1/  / 


vV/  =  v*-Vz-v^-Vi-i 


m-l  ,  .  . 

^j+m  ~  %+m-l  +  (— l)*^V*’/j-(.Tn-l  f  dt' ( 

k=0  .  \  '^ 

tj+m-l 

The  integrals  were  evaluated  euid  the  result  was  the  recursive  relation 

11  1 
7fc  +  +  37fc-2  H - h  ^TjTY'T'o  =  1 

for  fc  =  0, 1,  •  •  • ,  m  —  1  and  the  scheme  emerges  in  the  form 

m-l 

‘^j+m  =  Uj+fn-1  +  h  'yk^^ fj+m—l 

fc=0 


(18) 


(19) 
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Thet^cationerrorforthisclassofexplidtschemeswasgivenbyOrA”*)  Th,  ■  , 

m  -  2  was  the  example  important  for  the  solution  of  tlS  M  cl^,  ‘ 

scheme  was  then  regarded  either  as  an  exnlirit  <=  >,  ^  Navier-Stokes  equations.  This 

The  coefflcieme  j,  were  LlT,  -  >  Id  i  T  ”  “‘'P  =>  PC-scheme. 

-I  >  71  -  2  ana  the  scheme  was  given  by 

«,+2  =  +  ^h(3fi+i  -  fj)  +  0(e)  J20) 

It  requires  three  levels  of  the  solution  u,r.,  and  u  Th  u 

equation  The  scheme  applied  to  the 

du  du 

can  be  shown  (using  the  standard  Fourier  anaivsut  t, 

and  Moin,  1991),  but  the  instabiUty  sets  in  so  slowly ‘th^'t“‘rM" 

Marns-Moulton  multi-stev  imnlir.it  sche.mps 

vairrid  /  .ote‘“  r 

3+m  were  added  to  the  interpolation  polynomial 


with  s  defined  now  by 


leading  to 


s  = 


t'-t 


'j^Tn 


h 


m 

«i+m  =  tij+m-1  +  ^ 

Jfc=0 


(21) 

:z,'3arkisz“K  "rzr  “• 

1  1  1  1 

~I”  n^k — 1  ”f”  _ 2  “1“  *  *  "  "f"  _ ~1Jr\  n 

2  3*21-  ^22) 

for  A:  =  0,  •  • . ,  m.  The  scheme  was  implicit  since  f  -^  -  fd  \  j  . 

error  was  given  by  0(h’”+i)  The  sneHal  />oo  •' tij+mj  and  its  truncation 

Navier-Stote  equations.  vJ.  om3  I  wTan”  I  l"!l 

2 

«i+i  =  %  +  +  f.)  +  o(h^) 

It  required  iteration  to  compute  the  solution  at  the  new  time  level. 
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Multi-step  PC  (predictor-corrector)  schemes. 

The  explicit  Adams-Bashforth  and  the  implicit  Adams-Moulton  multi-step  schemes  can 
be  combined  into  a  system  of  predictor-corrector  (PC)  schemes.  The  general  form  of  the 
PC  schemes  followed  from  (19)  and  (21) 


m— 1 

+  h  +  0{h'^)  (23) 

Jb=0 

m 

'^i+m  —  1  “f*  h  ^  ^  (24) 

fc=0 

If  the  corrector  step  was  adjusted  to  the  same  truncation  error  as  the  predictor  step  the 
scheme 

m— 1 

=  Uj  +  hJ2  7fcV''/(tj,  Uj)  (23) 

k=0 

m— 1 

Uj+i  =Uj  +  hJ2  V*=/(tj+i,  (24) 

fc=0 

of  order  0(/i"*)  emerged.  The  special  case  m  =  2  with  truncation  error  0{h'^)  is  of 
particular  importance  and  leads  to  the  predictor  step 

^i+i  ~  "I"  ~  /(^i-i)%-i)]  (25) 

and  the  corrector  step 

Uj+i  =  Uj  +  ^h[f{tj+i,u*j^i)  +  f{tj,  Uj)]  (26) 

This  scheme  can  be  applied  to  the  Navier-Stokes  equations  if  the  right  side  f{t,u)  was 
interpreted  as  the  svun  of  convective,  pressure  gradient  and  viscous  terms.  The  stability 
properties  of  the  PC  scheme  (25),  (26)  were  superior  to  the  explicit  Adams-Bashforth 
scheme  (Grigorieff  vol.2,  1977,  section  2.10.4).  The  main  disadvantages  of  multi-step 
schemes  were  the  necessity  to  design  a  single  step  start-up  scheme  and  the  fact  that  the 
schemes  discussed  here  were  only  valid  for  constant  spa.cing  in  time. 
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Multi-Step  time  integration  schemes  for  the  Navier-Stokes  equations 

The  application  of  time  integration  schemes  developed  for  systems  of  odes  to  the  Navier- 
Stokes  equations  was  not  trivial  since  the  incompressible  equations  (1)  to  (3)  contained  the 
mass  balance  which  was  not  an  evolution  equation.  The  scheme  (25),  (26)  was  adapted  to 
the  Navier-Stokes  equations.  Denoting  the  spatial  discretizations  with  and  which 
will  be  specified  later,  the  momentum  balances  (2)  and  (3)  were  given  by 


du 

dt 


fx{u,v,p,Re) 


dv 

dt 


fy{u,v,p,Re) 


where 


and 


f  /  r>  X  5u  5u  Sp  1  ,S^u  S^u 

fx{u,v,p,Re)  =  -u- - V- - ^  -f  — (-jr^  -t-  -— -) 

Sx  5y  Sx  Re^Sx^  Sy^ 


fy{u,v,p,Re)  =  —u 


5v  5v 


5x  5y 


^  1  S'^v 

Sy^ 


The  predictor  step  consisted  then  of  two  operations,  first  the  velocity  is  predicted  according 
to  (25,  superscript  refers  to  time  level) 


+  ^At[3f^(u^,v'^,p'^,Re)  -  fx{u^-\v^-\p^-\Re)] 

and 

^*n+i  =  -y”  +  ^At[Zfy{u^,v^,p^,Re)  -  fy{u^-\v^-\p'^-\Re)] 

and  second,  the  predicted  velocity  was  then  updated  using  mass  balance.  The  velocity 
components  did  not  satisfy  the  discretized  mass  balance 


Sx  5y 


=  £)*"+!  ^  0 


The  updated  velocity  components  were  required  to  satisfy  mass  balance,  which  is  accom¬ 
plished  by 


Sx 

=  ^*”+1  _ 

Sy 

and 

Su'^+^ 

Sx  ^  Sy-  ^ 

leading  to 

S^^  <52$  , 

Sx^  '^S^~ 
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The  updated  velocity  components  must  satisfy  the  momentum  balances,  which 

leads  to  the  pressure  correction 


p/n+l 

^  ^  2At 


This  concluded  the  predictor  step.  It  was  important  to  note  at  this  point,  that  the  pressure 
correction  method  as  described  above  only  satisfies  one  of  the  two  boundary  conditions 
at  fixed  walls,  namely  the  condition  of  no  flow  normal  to  the  wall.  The  remaining  no-slip 
condition  can  be  enforced  in  the  predictor  step,  where  the  momentum  balances  were  solved 
for  the  primed  velocity  components.  It  must  be  tested  whether  this  approach  was  suitable 
for  highly  imsteady  flows  around  airfoils.  Both  no-slip  conditions  can  be  enforced  in  the 
pressure  correction  method  at  the  expense  of  an  additional  Poisson  equation. 

The  corrector  step  involves  also  two  operations  similar  to  the  predictor  step.  The 
updated  velocity  and  pressure  enter  the  predictor  step  according  to  (26)  as  follows 

=  w”  -b  -f  U{u^,v^,p^,Re)] 

^"n+i  ^  ^  -b 

The  velocity  components  and  the  pressure  must  be  updated  as  for  the 

predictor  step  to  enforce  the  mass  balance.  Mass  balance  required 

5x  ^  5y  ^ 

Hence,  the  correction  was  necessary 

Sx 


V 


=  ■0””+^ _ 

5y 


to  insmre  it,  which  leads  to 


SH*  .  .. 

+  =  -D 


Sx^  Sy^ 

where  denotes  the  divergence  of  The  pressure  associated  with  tlic 

updated  velocity  field  can  be  determined  as  follows.  Updated  momentum  balances  givini 
by 

=U^  +  ^At[U{u'^+\v'^+\p'^+^  +Pc,Re)  +  /.(n",u”,p-,He)] 


u 


and 


^n+l  =  ^  yn+1  + 


were  used  to  compute  the  pressure  correction Pc-  Subtracting  the  corrector  step  moriK'iit  uni 
balance  from  the  updated  versions  led  to 


Pc  = 


2$* 

At 
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and  the  updated  pressure  was  then  given  by 


It  was  clear  that  both  corrector  and  predictor  steps  required  updating  to  satisfy  mass 
balance.  The  question  remains  how  the  stability  of  this  type  of  scheme  is  modified  by  the 
updating  steps. 

Runge-Kutta  explicit  schemes 

The  family  of  Runge-Kutta  schemes  can  be  applied  to  the  unsteady  Navier-Stokes  equa¬ 
tions,  if  proper  care  was  taken  of  mass  balance.  The  momentum  balances  were  written  as 
in  the  previous  section  as 

du 

—  =  U{u,v,p,Re) 
dv 

—  =  fy{u,v,p,Re) 

together  with  mass  balance 

Ux+Vy=0 

The  standard  fourth  order  accmate  Runge-Kutta  method  was  considered  for  this  system 
of  equations.  The  solution  at  the  new  time  level  level  was  constructed  as  the  hnear  combi¬ 
nation  of  four  increments.  Mass  balance  must  be  satisfied  at  the  new  time  level,  which  can 
be  insured  by  requiring  it  for  each  of  the  increments.  The  increments  were  then  computed 
as  follows. 

Step  (1):  The  increments  fci(w)  and  kl(v)  were  computed  using  the  momentum  bal¬ 
ances  and  then  modified  to  enforce  mass  balance.  The  momentum  balances  produce  the 
increments  fci(tt)  and  ki(v) 


ki(u)  = 

ki(v)  =  Atfy{u^,v'^,p^,Re) 
which  did  not  satisfy  mass  balance 

-I-  ^kiiv)  =  #  0 

The  correction  potential  was  required  to  correct  the  velocity  increments 

ki{u)  =  ki(u)  - 

ki(v)  =  ki(v)  - 

dp 

It  was  determined  by  enforcing  mass  balance  for  the  corrected  increments,  which  led  to 


The  pressure  correction  was  then  obtained  from  the  requirement  that  the  corrected  incre¬ 
ments  satisfy  the  momentum  balances 

ki{u)  =  +p^^\Re) 

fcl(t.)  =  A(/„K,„”,p”+pW,iJe) 

leading  to 


The  increments  for  the  velocity  components  and  the  pressure  were  thus  computed. 

Step  (2):  The  increments  k2{u)  and  k2{v)  were  computed  using  the  momentum  bal¬ 
ances  and  then  modified  as  in  step  (1).  The  increments 

k(u)  =  Atf^(u^  +  ^ki(u),v'*  +  ^ki(v),p’^  -1-  ^pi^\Be) 

h(v)  =  ^tfy{u^  +  ifci («),!;"  +  -f-  ^pl}\Re) 

did  not  satisfy  meiss  balance 

^h{u)  +  ^k{v)  =  #  0 

and  the  correction  potential  was  required  as  before 

k2{u)  =  k2{u)  — 

ox 

k2{v)  =  k2(v)  — 

oy 

It  was  determined  by 

'•Sx^  Sy^' 

The  pressure  correction  was  then  obtained  from  the  requirement  that  the  corrected  incre¬ 
ments  satisfy  the  momentmn  balances 


kiiu)  =  A(/.(u”  +  ifciW.v"  +  +pf  ,Be) 


k2(v)  =  A(/„(ii”+  i*:, («),!,”  +  +pf\Re) 


leading  to 


p$(2) 

At 


The  increments  for  the  velocity  components  and  the  pressure  were  thus  computed. 
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Step  (3):  The  increments  k3{u)  and  k3{v)  were  computed  using  the  momentum  bal- 
ances  and  then  modified  as  in  step  (2).  The  increments 

kiu)  =  AtMu’'  +  ^k2{u):v”  +  +  \p?\Re) 

kz{v)  =  +  -k2{u),v^  +  ^k2{v),p^  +  ~p^^\Re) 

did  not  satisfy  mass  balance 


h{u)  +  -^hiv)  =  0 

oy 


Sx 


and  the  correction  potential  was  required  as  before 


kziu)  =  k{u)  -  ^$(3) 


It  was  determined  by 


k3(v)  =  k3{v)  — 

oy 


The  pressure  correction  was  then  obtained  from  the  requirement  that  the  corrected  incre 
ments  satisfied  the  momentum  balances 


hiu)  =  Atf^{u^  +  iA:2(u),u”  +  +p^^\Re) 


leading  to 


kziv)  =  Atfy{vr  +  \k2{u),v^  +  ]:k2{v),p^  +p(^\Re) 


p$(3) 

At 


The  increments  for  the  velocity  components  and  the  pressure  were  thus  computed. 

Step  (4).  The  increments  k^{u)  and  were  computed  using  the  momentum  bal¬ 
ances  and  then  modified  as  in  step  (3).  The  increments 


k^{u)  =  At/x(u"  -H  fc3(M),  v”  +  k3{v),p^  +p(^)^Re) 

hiv)  =  Atfy{u^  +  k3{u),  u”  -I-  k3{v),p^  +pf\Re) 
did  not  satisfy  mass  balance 


s  s 

k4{u)  4-  —k^{v)  =  ^  0 


5x 
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and  the  correction  potential  was  required  as  before 


k^{u)  =  k^{u)  - 
k^{v)  =  k^{v)  — 

oy 

It  was  determined  by 

The  pressure  correction  was  then  obtained  from  the  requirement  that  the  corrected  incre¬ 
ments  satisfy  the  momentum  balances 


k4iu)  =  AtU{u^  -f  k3{u),v^  +  h{v),p^  -h p^^\Re) 


leading  to 


k^iv)  =  Atfy{u'^  +ki{u),v'^  +  k3{v),p'^  +p^^\  Re) 


At 


The  increment  computation  was  thus  completed. 

The  velocity  and  pressure  at  the  new  time  level  were  then  given  by 


—  «”  +  -{ki{u)  +  2k2{u)  4-  2k3(u)  -1-  k4{u)) 
yn+l  ^  ^  1  ^  2A:2(v)  +  2k3iv)  +  k^iv)) 


=P^  +  i(p(i)  -1-  2p(2)  + 


according  to  the  standard  Runge-Kutta  method.  The  new  velocity  field  clearly  satisfies 
mass  balance  since  all  increments  do.  The  new  pressure  could  be  computed  as  = 
p”  +pi^^  or  as  linear  combination  of  the  pressure  corrections  associated  with  the  individual 
velocity  increments.  The  forth  order  Runge-Kutta  method  allows  easy  control  of  the  time 
step  according  to 


K  ^gp3(^)-fc2(n)|| 


<  K, 


max 


and 


<2 


Hfc3(^)  -fe2(^)i 
p2('y)  -  A:i(v)j 


<  K 


max 


where  the  maximum  norm  was  used  for  the  velocity  increments.  The  strategy  for  the  time 
step  control  was  then  to  double  the  time  step  ifrboth  ratios  were  below  Kmin  =  0.01  and 
to  half  the  time  step  if  one  of  the  two  ratios  was  above  =  0.05. 
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VORTICITY  MAGNITUDE 


1 


0.0:  Vorticity  magnitude  at  time  t  =  898. 


VORTICITY  MAGNITUDE 


113  uM  vap.  Pentane,  60’  Heated  Jet,  Ti=305,  Cs=.05,  Sch=l. 

300-1 - ^ ^ ^ ^ - , 


[^^uim]  uoisjadsiQ 
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Figure  3:  Spatial  dispersion  of  113  fxm  vaporizing  Pentane  droplets  in  a  heated  turbulent 
round  jet  (entrance  temperature  Ti  =  305"iir,  LES  simulation  with  Cs  =  0.05,  Sc  =  1.0). 


HexaDecane  Particles, Tj=20’,x/D=30,ils=9ms 
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Figure  4:  Axial  velocity  auto-correlation  function  for  50  //m  Hexadecane  droplets. 


60  D  LES,  113  uM  vaporizing  pentane  particles 
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Figure  7:  Dispersion  of  water  droplets  in  a  particle  laden  turbulent  round  jet  for  different 
mass  loadings  (experiments:  left  graph)  and  the  mass  loading  5%  (LES  simulation:  right 
graph).  Top:  Photo  of  Spray. 


Burger’s  Equation,  Re  =  100. 
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Figure  13:  Time  convergence  of  the  third  order  Rimge-Kutta  time  integration  scheme 
and  6th  order  compact  spatial  discretization  for  the  Burgers  equation  {Re  =  100). 


Convergence  with  grid  refinement 
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In(error) 


Burger’s  Equation,  Re  =  100. 


Convergence  with  arid  refinement 
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Figure  14:  Spatial  convergence  of  the  third  order  Runge-Kutta  time  integration  sc! 
and  6th  order  compact  spatial  discretization  for  the  Burgers  equation  {Re  =  100). 


VORTiaTY  MAGNITUDE 


Figure  15:  Vorticity  magnitude  for  an  axisymmetric  round  jet  using  the  third  order 
Runge-Kutta  time  integration  scheme  at  Re  =  5000  and  t  =  46. 
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Figure  18:  LES  simulation  of  a  plane  jet  using  the  fourth  order  Runge-Kutta  time  integra¬ 
tion  scheme  and  fourth  order  compact  spatial  discretization,  velocity  vectors  for  Re  =  2500 
at  t  =  25. 


I:  LES  simulation  of  the  backward-facing  step  flow  using  the  fourth  order  Runge- 
»!  integration  scheme,  iso-lines  of  vorticity  magnitude  at  t  =  15  and  Re  =  500( 
step  height). 


Figure  21:  Spray  size  distribution. 
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Figure  22:  Apparatus  for  tracking  fluorescent  particles  in  a  turbulent  spray. 


